plus sc10x

run cellranger with Spp1-EGFP added reference index

load dependancies

Spp1-GFP mice, SIMPLE induced stroke,
P28 7d: GFP(Spp1+) x2, MIG(Microglia) x2, CTR(Contra) x1
P05 3d: GFP(Spp1+) x2, MIG(Microglia) x2, CTR(Contra) x1

loading 70k cells, CX3CR1+
demo, cellranger called 44,614 cells
plus with Spp1-EGFP, cellranger called 47,886 cells

load 10x data

filt.10x <- Read10X(data.dir = "../output_plus_exogene/filtered_feature_bc_matrix/")
## 10X data contains more than one type and is being returned as a list containing matrices of each type.
GEX <- filt.10x$`Gene Expression`
FB <- filt.10x$`Antibody Capture`

check datasets

dim(GEX)
## [1] 32286 47886
GEX[1:6,1:6]
## 6 x 6 sparse Matrix of class "dgCMatrix"
##         AAACCCAAGAAATCCA-1 AAACCCAAGAATCGTA-1 AAACCCAAGAGCCTGA-1
## Xkr4                     .                  .                  .
## Gm1992                   .                  .                  .
## Gm19938                  .                  .                  .
## Gm37381                  .                  .                  .
## Rp1                      .                  .                  .
## Sox17                    .                  .                  .
##         AAACCCAAGATGCTAA-1 AAACCCAAGCACCTGC-1 AAACCCAAGCGATCGA-1
## Xkr4                     .                  .                  .
## Gm1992                   .                  .                  .
## Gm19938                  .                  .                  .
## Gm37381                  .                  .                  .
## Rp1                      .                  .                  .
## Sox17                    .                  .                  .
dim(FB)
## [1]    10 47886
FB[,1:6]
## 10 x 6 sparse Matrix of class "dgCMatrix"
##           AAACCCAAGAAATCCA-1 AAACCCAAGAATCGTA-1 AAACCCAAGAGCCTGA-1
## P28.GFP.1                 12                 49                392
## P28.GFP.2                 26                 10                 29
## P28.MIG.1                 11                  9                 21
## P28.MIG.2                  5                  1                  1
## P28.CTR.1                  4                  .                  1
## P05.GFP.1                139                  .                  4
## P05.GFP.2                 15                  3                  7
## P05.MIG.1                  2                  4                  2
## P05.MIG.2                  2                  2                209
## P05.CTR.1                  5                  2                  4
##           AAACCCAAGATGCTAA-1 AAACCCAAGCACCTGC-1 AAACCCAAGCGATCGA-1
## P28.GFP.1                  3                 23                 10
## P28.GFP.2                 28                 68                  6
## P28.MIG.1                  6                 23                  9
## P28.MIG.2                 51                  7                  9
## P28.CTR.1                  .                  5                  .
## P05.GFP.1                  4                 10                  3
## P05.GFP.2                  3                266                  4
## P05.MIG.1                  3                 78                  4
## P05.MIG.2                  4                259                188
## P05.CTR.1                 40                  8                  4
rowSums(FB)
## P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1 
##   1648415   4014042   1456711    974485    312662   1004996   1222480    876714 
## P05.MIG.2 P05.CTR.1 
##   1141690    824693
rowSums(FB>0)
## P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1 
##     47701     47859     47690     46474     35812     46092     47378     44721 
## P05.MIG.2 P05.CTR.1 
##     46311     42039

FB

# build seurat object
FB.seur <- CreateSeuratObject(counts = FB,project = "plus_SIMPLE")

FB.seur
## An object of class Seurat 
## 10 features across 47886 samples within 1 assay 
## Active assay: RNA (10 features, 0 variable features)
rownames(FB.seur)
##  [1] "P28.GFP.1" "P28.GFP.2" "P28.MIG.1" "P28.MIG.2" "P28.CTR.1" "P05.GFP.1"
##  [7] "P05.GFP.2" "P05.MIG.1" "P05.MIG.2" "P05.CTR.1"

normalization

perform Seurat ’Demultiplexing with hashtag oligos(HTOs)
ref https://satijalab.org/seurat/articles/hashing_vignette.html

# normalize
FB.seur <- NormalizeData(FB.seur)
FB.seur <- FindVariableFeatures(FB.seur, selection.method = "mean.var.plot")
FB.seur <- ScaleData(FB.seur, features = VariableFeatures(FB.seur))
## Centering and scaling data matrix
# independent assay for HTO data  
FB.seur[['HTO']] <- CreateAssayObject(counts = FB)
FB.seur <- NormalizeData(FB.seur, assay = 'HTO', normalization.method = 'CLR')
## Normalizing across features

check demux

first define FB colors based on conditions

scales::show_col(ggsci::pal_igv("default")(49))

color.FB <- ggsci::pal_igv("default")(49)[c(2,13,33,1,15,
                                            34,26,28,40,18)]


level.FB <- c(rownames(FB.seur))

color.FBraw <-  c("#FF6C91","lightgrey",color.FB)
scales::show_col(color.FBraw, ncol = 5)

scales::show_col(color.FB, ncol = 5)

par(mfrow=c(2,5))

#plot(sort(rowSums(t(FB))),
#     xlab="reordered cells", 
#     ylab="UMI counts", log="xy",
#     main="sum")


plot(sort(t(FB)[,1]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[1])
plot(sort(t(FB)[,2]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[2])
plot(sort(t(FB)[,3]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[3])
plot(sort(t(FB)[,4]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[4])
plot(sort(t(FB)[,5]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[5])
plot(sort(t(FB)[,6]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[6])
plot(sort(t(FB)[,7]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[7])
plot(sort(t(FB)[,8]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[8])
plot(sort(t(FB)[,9]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[9])
plot(sort(t(FB)[,10]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[10])

range ‘q’

q0.50 ~ 0.99
FB.test <- FB.seur
#cut.list <- list()

table.demux <- data.frame(sample=c("quantile","Doublet",rownames(FB.test),"Negative"))

for(i in 50:99){
  FB.test <- HTODemux(FB.test, assay = "HTO", positive.quantile = i/100)
  table.demux <- cbind(table.demux,
                       c(i/100,table(FB.test$hash.ID)[c("Doublet",rownames(FB.test),"Negative")]))
}
## Cutoff for P28.GFP.1 : 7 reads
## Cutoff for P28.GFP.2 : 15 reads
## Cutoff for P28.MIG.1 : 7 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 7 reads
## Cutoff for P05.MIG.1 : 2 reads
## Cutoff for P05.MIG.2 : 4 reads
## Cutoff for P05.CTR.1 : 2 reads
## Cutoff for P28.GFP.1 : 7 reads
## Cutoff for P28.GFP.2 : 15 reads
## Cutoff for P28.MIG.1 : 8 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 7 reads
## Cutoff for P05.MIG.1 : 2 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 2 reads
## Cutoff for P28.GFP.1 : 7 reads
## Cutoff for P28.GFP.2 : 16 reads
## Cutoff for P28.MIG.1 : 8 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 7 reads
## Cutoff for P05.MIG.1 : 2 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 2 reads
## Cutoff for P28.GFP.1 : 7 reads
## Cutoff for P28.GFP.2 : 16 reads
## Cutoff for P28.MIG.1 : 8 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 7 reads
## Cutoff for P05.MIG.1 : 2 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 2 reads
## Cutoff for P28.GFP.1 : 7 reads
## Cutoff for P28.GFP.2 : 16 reads
## Cutoff for P28.MIG.1 : 8 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 7 reads
## Cutoff for P05.MIG.1 : 2 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 8 reads
## Cutoff for P28.GFP.2 : 16 reads
## Cutoff for P28.MIG.1 : 8 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 8 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 8 reads
## Cutoff for P28.GFP.2 : 16 reads
## Cutoff for P28.MIG.1 : 9 reads
## Cutoff for P28.MIG.2 : 4 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 8 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 8 reads
## Cutoff for P28.GFP.2 : 17 reads
## Cutoff for P28.MIG.1 : 9 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 8 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 5 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 8 reads
## Cutoff for P28.GFP.2 : 17 reads
## Cutoff for P28.MIG.1 : 9 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 3 reads
## Cutoff for P05.GFP.2 : 8 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 8 reads
## Cutoff for P28.GFP.2 : 17 reads
## Cutoff for P28.MIG.1 : 9 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 8 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 9 reads
## Cutoff for P28.GFP.2 : 17 reads
## Cutoff for P28.MIG.1 : 10 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 8 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 9 reads
## Cutoff for P28.GFP.2 : 18 reads
## Cutoff for P28.MIG.1 : 10 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 9 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 9 reads
## Cutoff for P28.GFP.2 : 18 reads
## Cutoff for P28.MIG.1 : 10 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 9 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 9 reads
## Cutoff for P28.GFP.2 : 18 reads
## Cutoff for P28.MIG.1 : 10 reads
## Cutoff for P28.MIG.2 : 5 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 9 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 9 reads
## Cutoff for P28.GFP.2 : 18 reads
## Cutoff for P28.MIG.1 : 10 reads
## Cutoff for P28.MIG.2 : 6 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 9 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 6 reads
## Cutoff for P05.CTR.1 : 3 reads
## Cutoff for P28.GFP.1 : 10 reads
## Cutoff for P28.GFP.2 : 18 reads
## Cutoff for P28.MIG.1 : 11 reads
## Cutoff for P28.MIG.2 : 6 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 9 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 7 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 10 reads
## Cutoff for P28.GFP.2 : 19 reads
## Cutoff for P28.MIG.1 : 11 reads
## Cutoff for P28.MIG.2 : 6 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 10 reads
## Cutoff for P05.MIG.1 : 3 reads
## Cutoff for P05.MIG.2 : 7 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 10 reads
## Cutoff for P28.GFP.2 : 19 reads
## Cutoff for P28.MIG.1 : 11 reads
## Cutoff for P28.MIG.2 : 6 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 10 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 7 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 10 reads
## Cutoff for P28.GFP.2 : 19 reads
## Cutoff for P28.MIG.1 : 12 reads
## Cutoff for P28.MIG.2 : 6 reads
## Cutoff for P28.CTR.1 : 1 reads
## Cutoff for P05.GFP.1 : 4 reads
## Cutoff for P05.GFP.2 : 10 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 7 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 10 reads
## Cutoff for P28.GFP.2 : 19 reads
## Cutoff for P28.MIG.1 : 12 reads
## Cutoff for P28.MIG.2 : 6 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 10 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 7 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 11 reads
## Cutoff for P28.GFP.2 : 20 reads
## Cutoff for P28.MIG.1 : 12 reads
## Cutoff for P28.MIG.2 : 7 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 11 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 8 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 11 reads
## Cutoff for P28.GFP.2 : 20 reads
## Cutoff for P28.MIG.1 : 12 reads
## Cutoff for P28.MIG.2 : 7 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 11 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 8 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 11 reads
## Cutoff for P28.GFP.2 : 20 reads
## Cutoff for P28.MIG.1 : 13 reads
## Cutoff for P28.MIG.2 : 7 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 11 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 8 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 11 reads
## Cutoff for P28.GFP.2 : 21 reads
## Cutoff for P28.MIG.1 : 13 reads
## Cutoff for P28.MIG.2 : 7 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 11 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 8 reads
## Cutoff for P05.CTR.1 : 4 reads
## Cutoff for P28.GFP.1 : 12 reads
## Cutoff for P28.GFP.2 : 21 reads
## Cutoff for P28.MIG.1 : 13 reads
## Cutoff for P28.MIG.2 : 7 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 11 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 8 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 12 reads
## Cutoff for P28.GFP.2 : 21 reads
## Cutoff for P28.MIG.1 : 14 reads
## Cutoff for P28.MIG.2 : 7 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 12 reads
## Cutoff for P05.MIG.1 : 4 reads
## Cutoff for P05.MIG.2 : 9 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 12 reads
## Cutoff for P28.GFP.2 : 21 reads
## Cutoff for P28.MIG.1 : 14 reads
## Cutoff for P28.MIG.2 : 8 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 5 reads
## Cutoff for P05.GFP.2 : 12 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 9 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 13 reads
## Cutoff for P28.GFP.2 : 22 reads
## Cutoff for P28.MIG.1 : 14 reads
## Cutoff for P28.MIG.2 : 8 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 6 reads
## Cutoff for P05.GFP.2 : 12 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 9 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 13 reads
## Cutoff for P28.GFP.2 : 22 reads
## Cutoff for P28.MIG.1 : 15 reads
## Cutoff for P28.MIG.2 : 8 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 6 reads
## Cutoff for P05.GFP.2 : 13 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 9 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 13 reads
## Cutoff for P28.GFP.2 : 22 reads
## Cutoff for P28.MIG.1 : 15 reads
## Cutoff for P28.MIG.2 : 8 reads
## Cutoff for P28.CTR.1 : 2 reads
## Cutoff for P05.GFP.1 : 6 reads
## Cutoff for P05.GFP.2 : 13 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 10 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 14 reads
## Cutoff for P28.GFP.2 : 23 reads
## Cutoff for P28.MIG.1 : 16 reads
## Cutoff for P28.MIG.2 : 9 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 6 reads
## Cutoff for P05.GFP.2 : 13 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 10 reads
## Cutoff for P05.CTR.1 : 5 reads
## Cutoff for P28.GFP.1 : 14 reads
## Cutoff for P28.GFP.2 : 23 reads
## Cutoff for P28.MIG.1 : 16 reads
## Cutoff for P28.MIG.2 : 9 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 6 reads
## Cutoff for P05.GFP.2 : 14 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 10 reads
## Cutoff for P05.CTR.1 : 6 reads
## Cutoff for P28.GFP.1 : 14 reads
## Cutoff for P28.GFP.2 : 23 reads
## Cutoff for P28.MIG.1 : 17 reads
## Cutoff for P28.MIG.2 : 9 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 6 reads
## Cutoff for P05.GFP.2 : 14 reads
## Cutoff for P05.MIG.1 : 5 reads
## Cutoff for P05.MIG.2 : 10 reads
## Cutoff for P05.CTR.1 : 6 reads
## Cutoff for P28.GFP.1 : 15 reads
## Cutoff for P28.GFP.2 : 24 reads
## Cutoff for P28.MIG.1 : 17 reads
## Cutoff for P28.MIG.2 : 9 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 7 reads
## Cutoff for P05.GFP.2 : 14 reads
## Cutoff for P05.MIG.1 : 6 reads
## Cutoff for P05.MIG.2 : 11 reads
## Cutoff for P05.CTR.1 : 6 reads
## Cutoff for P28.GFP.1 : 15 reads
## Cutoff for P28.GFP.2 : 24 reads
## Cutoff for P28.MIG.1 : 18 reads
## Cutoff for P28.MIG.2 : 10 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 7 reads
## Cutoff for P05.GFP.2 : 15 reads
## Cutoff for P05.MIG.1 : 6 reads
## Cutoff for P05.MIG.2 : 11 reads
## Cutoff for P05.CTR.1 : 6 reads
## Cutoff for P28.GFP.1 : 16 reads
## Cutoff for P28.GFP.2 : 25 reads
## Cutoff for P28.MIG.1 : 18 reads
## Cutoff for P28.MIG.2 : 10 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 7 reads
## Cutoff for P05.GFP.2 : 15 reads
## Cutoff for P05.MIG.1 : 6 reads
## Cutoff for P05.MIG.2 : 11 reads
## Cutoff for P05.CTR.1 : 6 reads
## Cutoff for P28.GFP.1 : 16 reads
## Cutoff for P28.GFP.2 : 25 reads
## Cutoff for P28.MIG.1 : 19 reads
## Cutoff for P28.MIG.2 : 11 reads
## Cutoff for P28.CTR.1 : 3 reads
## Cutoff for P05.GFP.1 : 7 reads
## Cutoff for P05.GFP.2 : 16 reads
## Cutoff for P05.MIG.1 : 6 reads
## Cutoff for P05.MIG.2 : 12 reads
## Cutoff for P05.CTR.1 : 7 reads
## Cutoff for P28.GFP.1 : 16 reads
## Cutoff for P28.GFP.2 : 26 reads
## Cutoff for P28.MIG.1 : 19 reads
## Cutoff for P28.MIG.2 : 11 reads
## Cutoff for P28.CTR.1 : 4 reads
## Cutoff for P05.GFP.1 : 7 reads
## Cutoff for P05.GFP.2 : 16 reads
## Cutoff for P05.MIG.1 : 6 reads
## Cutoff for P05.MIG.2 : 12 reads
## Cutoff for P05.CTR.1 : 7 reads
## Cutoff for P28.GFP.1 : 17 reads
## Cutoff for P28.GFP.2 : 26 reads
## Cutoff for P28.MIG.1 : 20 reads
## Cutoff for P28.MIG.2 : 11 reads
## Cutoff for P28.CTR.1 : 4 reads
## Cutoff for P05.GFP.1 : 8 reads
## Cutoff for P05.GFP.2 : 17 reads
## Cutoff for P05.MIG.1 : 7 reads
## Cutoff for P05.MIG.2 : 13 reads
## Cutoff for P05.CTR.1 : 7 reads
## Cutoff for P28.GFP.1 : 18 reads
## Cutoff for P28.GFP.2 : 27 reads
## Cutoff for P28.MIG.1 : 21 reads
## Cutoff for P28.MIG.2 : 12 reads
## Cutoff for P28.CTR.1 : 4 reads
## Cutoff for P05.GFP.1 : 8 reads
## Cutoff for P05.GFP.2 : 17 reads
## Cutoff for P05.MIG.1 : 7 reads
## Cutoff for P05.MIG.2 : 13 reads
## Cutoff for P05.CTR.1 : 7 reads
## Cutoff for P28.GFP.1 : 18 reads
## Cutoff for P28.GFP.2 : 27 reads
## Cutoff for P28.MIG.1 : 22 reads
## Cutoff for P28.MIG.2 : 12 reads
## Cutoff for P28.CTR.1 : 4 reads
## Cutoff for P05.GFP.1 : 8 reads
## Cutoff for P05.GFP.2 : 18 reads
## Cutoff for P05.MIG.1 : 7 reads
## Cutoff for P05.MIG.2 : 14 reads
## Cutoff for P05.CTR.1 : 8 reads
## Cutoff for P28.GFP.1 : 19 reads
## Cutoff for P28.GFP.2 : 28 reads
## Cutoff for P28.MIG.1 : 23 reads
## Cutoff for P28.MIG.2 : 13 reads
## Cutoff for P28.CTR.1 : 5 reads
## Cutoff for P05.GFP.1 : 9 reads
## Cutoff for P05.GFP.2 : 18 reads
## Cutoff for P05.MIG.1 : 7 reads
## Cutoff for P05.MIG.2 : 14 reads
## Cutoff for P05.CTR.1 : 8 reads
## Cutoff for P28.GFP.1 : 20 reads
## Cutoff for P28.GFP.2 : 29 reads
## Cutoff for P28.MIG.1 : 24 reads
## Cutoff for P28.MIG.2 : 13 reads
## Cutoff for P28.CTR.1 : 5 reads
## Cutoff for P05.GFP.1 : 9 reads
## Cutoff for P05.GFP.2 : 19 reads
## Cutoff for P05.MIG.1 : 8 reads
## Cutoff for P05.MIG.2 : 15 reads
## Cutoff for P05.CTR.1 : 8 reads
## Cutoff for P28.GFP.1 : 21 reads
## Cutoff for P28.GFP.2 : 30 reads
## Cutoff for P28.MIG.1 : 25 reads
## Cutoff for P28.MIG.2 : 14 reads
## Cutoff for P28.CTR.1 : 5 reads
## Cutoff for P05.GFP.1 : 9 reads
## Cutoff for P05.GFP.2 : 20 reads
## Cutoff for P05.MIG.1 : 8 reads
## Cutoff for P05.MIG.2 : 16 reads
## Cutoff for P05.CTR.1 : 9 reads
## Cutoff for P28.GFP.1 : 22 reads
## Cutoff for P28.GFP.2 : 30 reads
## Cutoff for P28.MIG.1 : 26 reads
## Cutoff for P28.MIG.2 : 15 reads
## Cutoff for P28.CTR.1 : 6 reads
## Cutoff for P05.GFP.1 : 10 reads
## Cutoff for P05.GFP.2 : 21 reads
## Cutoff for P05.MIG.1 : 9 reads
## Cutoff for P05.MIG.2 : 17 reads
## Cutoff for P05.CTR.1 : 9 reads
## Cutoff for P28.GFP.1 : 23 reads
## Cutoff for P28.GFP.2 : 31 reads
## Cutoff for P28.MIG.1 : 27 reads
## Cutoff for P28.MIG.2 : 16 reads
## Cutoff for P28.CTR.1 : 6 reads
## Cutoff for P05.GFP.1 : 10 reads
## Cutoff for P05.GFP.2 : 22 reads
## Cutoff for P05.MIG.1 : 9 reads
## Cutoff for P05.MIG.2 : 18 reads
## Cutoff for P05.CTR.1 : 10 reads
## Cutoff for P28.GFP.1 : 24 reads
## Cutoff for P28.GFP.2 : 33 reads
## Cutoff for P28.MIG.1 : 29 reads
## Cutoff for P28.MIG.2 : 17 reads
## Cutoff for P28.CTR.1 : 7 reads
## Cutoff for P05.GFP.1 : 11 reads
## Cutoff for P05.GFP.2 : 23 reads
## Cutoff for P05.MIG.1 : 10 reads
## Cutoff for P05.MIG.2 : 19 reads
## Cutoff for P05.CTR.1 : 11 reads
## Cutoff for P28.GFP.1 : 26 reads
## Cutoff for P28.GFP.2 : 34 reads
## Cutoff for P28.MIG.1 : 32 reads
## Cutoff for P28.MIG.2 : 18 reads
## Cutoff for P28.CTR.1 : 7 reads
## Cutoff for P05.GFP.1 : 12 reads
## Cutoff for P05.GFP.2 : 25 reads
## Cutoff for P05.MIG.1 : 11 reads
## Cutoff for P05.MIG.2 : 20 reads
## Cutoff for P05.CTR.1 : 12 reads
## Cutoff for P28.GFP.1 : 29 reads
## Cutoff for P28.GFP.2 : 37 reads
## Cutoff for P28.MIG.1 : 35 reads
## Cutoff for P28.MIG.2 : 20 reads
## Cutoff for P28.CTR.1 : 8 reads
## Cutoff for P05.GFP.1 : 13 reads
## Cutoff for P05.GFP.2 : 27 reads
## Cutoff for P05.MIG.1 : 12 reads
## Cutoff for P05.MIG.2 : 22 reads
## Cutoff for P05.CTR.1 : 13 reads
## Cutoff for P28.GFP.1 : 33 reads
## Cutoff for P28.GFP.2 : 40 reads
## Cutoff for P28.MIG.1 : 41 reads
## Cutoff for P28.MIG.2 : 24 reads
## Cutoff for P28.CTR.1 : 10 reads
## Cutoff for P05.GFP.1 : 15 reads
## Cutoff for P05.GFP.2 : 31 reads
## Cutoff for P05.MIG.1 : 14 reads
## Cutoff for P05.MIG.2 : 26 reads
## Cutoff for P05.CTR.1 : 15 reads
rownames(table.demux) <- table.demux$sample
table.demux <- table.demux[,2:ncol(table.demux)]
table.demux <- data.frame(t(table.demux))
rownames(table.demux) <- NULL
table.demux
##    quantile Doublet P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1
## 1      0.50   46325       127       175       144       134       140       187
## 2      0.51   46090       156       205       152       167       158       219
## 3      0.52   45982       167       205       158       181       170       233
## 4      0.53   45982       167       205       158       181       170       233
## 5      0.54   45774       186       234       177       208       190       249
## 6      0.55   45113       216       314       244       285       272       337
## 7      0.56   45002       227       329       242       295       281       351
## 8      0.57   44553       276       356       285       312       323       415
## 9      0.58   44408       293       376       297       323       335       437
## 10     0.59   43900       333       441       354       377       391       433
## 11     0.60   43562       337       493       373       421       426       468
## 12     0.61   43266       374       497       397       450       471       495
## 13     0.62   43266       374       497       397       450       471       495
## 14     0.63   43266       374       497       397       450       471       495
## 15     0.64   43042       389       535       418       447       500       519
## 16     0.65   42430       435       610       459       507       568       596
## 17     0.66   42085       485       621       496       535       611       636
## 18     0.67   41454       543       702       552       597       688       721
## 19     0.68   41364       548       718       553       603       700       730
## 20     0.69   39706       689       912       739       809       766       857
## 21     0.70   38710       798       996       842       880       881       959
## 22     0.71   38710       798       996       842       880       881       959
## 23     0.72   38611       807      1009       840       888       893       971
## 24     0.73   38280       868      1008       879       915       927      1012
## 25     0.74   37899       886      1060       928       962       958      1055
## 26     0.75   37727       900      1082       940       977       976      1076
## 27     0.76   37023       965      1168      1021      1012      1053      1166
## 28     0.77   35964      1048      1260      1129      1124      1168      1220
## 29     0.78   35864      1055      1271      1133      1137      1175      1232
## 30     0.79   35806      1060      1277      1136      1144      1182      1237
## 31     0.80   34516      1180      1399      1262      1270      1249      1381
## 32     0.81   34318      1196      1423      1283      1295      1264      1402
## 33     0.82   34268      1199      1432      1282      1298      1272      1408
## 34     0.83   32832      1346      1581      1429      1445      1428      1516
## 35     0.84   32647      1357      1610      1441      1453      1447      1533
## 36     0.85   32198      1390      1629      1480      1501      1497      1585
## 37     0.86   31882      1406      1685      1508      1522      1540      1621
## 38     0.87   31215      1468      1732      1569      1596      1562      1705
## 39     0.88   30315      1531      1866      1655      1700      1666      1742
## 40     0.89   29804      1579      1898      1710      1735      1733      1791
## 41     0.90   29622      1589      1934      1724      1752      1744      1807
## 42     0.91   28647      1662      2030      1822      1855      1817      1854
## 43     0.92   27908      1728      2079      1908      1952      1892      1930
## 44     0.93   27372      1772      2111      1951      2001      1942      1990
## 45     0.94   26689      1833      2226      2019      2050      1990      2025
## 46     0.95   26146      1884      2250      2062      2089      2041      2094
## 47     0.96   24787      2011      2372      2207      2205      2141      2207
## 48     0.97   23929      2077      2460      2296      2281      2209      2282
## 49     0.98   22459      2215      2543      2427      2397      2315      2419
## 50     0.99   20831      2348      2655      2514      2512      2384      2570
##    P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 Negative
## 1        126       179       159       155       35
## 2        140       215       167       178       39
## 3        149       229       181       192       39
## 4        149       229       181       192       39
## 5        177       256       203       187       45
## 6        231       291       279       249       55
## 7        241       306       292       263       57
## 8        290       352       346       314       64
## 9        307       372       343       328       67
## 10       374       424       393       389       77
## 11       408       472       424       419       83
## 12       428       508       451       460       89
## 13       428       508       451       460       89
## 14       428       508       451       460       89
## 15       443       544       470       486       93
## 16       506       634       516       522      103
## 17       537       671       546       557      106
## 18       603       662       612       629      123
## 19       617       669       621       640      123
## 20       806       851       811       796      144
## 21       908       968       888       889      167
## 22       908       968       888       889      167
## 23       914       978       905       900      170
## 24       942      1017       932       931      175
## 25       985      1060       963       943      187
## 26       999      1082       974       954      199
## 27      1075      1092      1070      1020      221
## 28      1200      1233      1153      1146      241
## 29      1206      1253      1160      1156      244
## 30      1214      1265      1153      1160      252
## 31      1351      1427      1263      1302      286
## 32      1368      1456      1287      1294      300
## 33      1371      1465      1291      1298      302
## 34      1535      1551      1429      1456      338
## 35      1546      1575      1448      1481      348
## 36      1594      1626      1498      1531      357
## 37      1613      1655      1525      1550      379
## 38      1677      1728      1603      1624      407
## 39      1760      1785      1682      1720      464
## 40      1816      1818      1738      1779      485
## 41      1839      1841      1754      1777      503
## 42      1942      1949      1859      1879      570
## 43      2008      1990      1929      1959      603
## 44      2070      2050      1974      2008      645
## 45      2128      2102      2024      2082      718
## 46      2187      2163      2082      2120      768
## 47      2325      2300      2224      2226      881
## 48      2398      2378      2303      2309      964
## 49      2543      2526      2486      2446     1110
## 50      2683      2691      2626      2584     1488
#color.FB <- scales::hue_pal()(10)
par(mfrow=c(3,5))


plot(table.demux$Doublet, type="p", col="#FF6C91", pch=16, ylab="Doublet")
plot(table.demux$Negative, type="p", col="lightgrey", pch=16, ylab="Negative")
plot(rowSums(table.demux[,grep("P",colnames(table.demux))]), type="p", col="#306000", pch=16, ylab="Singlet")

plot(table.demux$quantile, pch=16, ylab="mod-quantile")
plot.new()

plot(table.demux[,2+1], type="p", col=color.FB[1], pch=16, ylab=colnames(table.demux)[2+1])
plot(table.demux[,2+2], type="p", col=color.FB[2], pch=16, ylab=colnames(table.demux)[2+2])
plot(table.demux[,2+3], type="p", col=color.FB[3], pch=16, ylab=colnames(table.demux)[2+3])
plot(table.demux[,2+4], type="p", col=color.FB[4], pch=16, ylab=colnames(table.demux)[2+4])
plot(table.demux[,2+5], type="p", col=color.FB[5], pch=16, ylab=colnames(table.demux)[2+5])

plot(table.demux[,2+6], type="p", col=color.FB[6], pch=16, ylab=colnames(table.demux)[2+6])
plot(table.demux[,2+7], type="p", col=color.FB[7], pch=16, ylab=colnames(table.demux)[2+7])
plot(table.demux[,2+8], type="p", col=color.FB[8], pch=16, ylab=colnames(table.demux)[2+8])
plot(table.demux[,2+9], type="p", col=color.FB[9], pch=16, ylab=colnames(table.demux)[2+9])
plot(table.demux[,2+10], type="p", col=color.FB[10], pch=16, ylab=colnames(table.demux)[2+10])

q0.9900 ~ 0.9999
#color.FB <- scales::hue_pal()(10)
par(mfrow=c(3,5))


plot(table.demux.s$Doublet, type="p", col="#FF6C91", pch=16, ylab="Doublet")
plot(table.demux.s$Negative, type="p", col="lightgrey", pch=16, ylab="Negative")
plot(rowSums(table.demux.s[,grep("P",colnames(table.demux.s))]), type="p", col="#306000", pch=16, ylab="Singlet")

plot(table.demux.s$quantile, pch=16, ylab="mod-quantile")
plot.new()

plot(table.demux.s[,2+1], type="p", col=color.FB[1], pch=16, ylab=colnames(table.demux.s)[2+1])
plot(table.demux.s[,2+2], type="p", col=color.FB[2], pch=16, ylab=colnames(table.demux.s)[2+2])
plot(table.demux.s[,2+3], type="p", col=color.FB[3], pch=16, ylab=colnames(table.demux.s)[2+3])
plot(table.demux.s[,2+4], type="p", col=color.FB[4], pch=16, ylab=colnames(table.demux.s)[2+4])
plot(table.demux.s[,2+5], type="p", col=color.FB[5], pch=16, ylab=colnames(table.demux.s)[2+5])

plot(table.demux.s[,2+6], type="p", col=color.FB[6], pch=16, ylab=colnames(table.demux.s)[2+6])
plot(table.demux.s[,2+7], type="p", col=color.FB[7], pch=16, ylab=colnames(table.demux.s)[2+7])
plot(table.demux.s[,2+8], type="p", col=color.FB[8], pch=16, ylab=colnames(table.demux.s)[2+8])
plot(table.demux.s[,2+9], type="p", col=color.FB[9], pch=16, ylab=colnames(table.demux.s)[2+9])
plot(table.demux.s[,2+10], type="p", col=color.FB[10], pch=16, ylab=colnames(table.demux.s)[2+10])

demultiplexing

(tags1/2/6-10 q0.998; tags3-5 q0.99)

# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.994)
## Cutoff for P28.GFP.1 : 36 reads
## Cutoff for P28.GFP.2 : 43 reads
## Cutoff for P28.MIG.1 : 45 reads
## Cutoff for P28.MIG.2 : 26 reads
## Cutoff for P28.CTR.1 : 12 reads
## Cutoff for P05.GFP.1 : 17 reads
## Cutoff for P05.GFP.2 : 34 reads
## Cutoff for P05.MIG.1 : 15 reads
## Cutoff for P05.MIG.2 : 29 reads
## Cutoff for P05.CTR.1 : 17 reads
table(FB.seur$HTO_classification.global)
## 
##  Doublet Negative  Singlet 
##    19627     1851    26408
table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]
## 
##   Doublet  Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 
##     19627      1851      2461      2741      2559      2613      2374      2672 
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 
##      2774      2782      2745      2687
# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.996)
## Cutoff for P28.GFP.1 : 38 reads
## Cutoff for P28.GFP.2 : 45 reads
## Cutoff for P28.MIG.1 : 48 reads
## Cutoff for P28.MIG.2 : 28 reads
## Cutoff for P28.CTR.1 : 13 reads
## Cutoff for P05.GFP.1 : 18 reads
## Cutoff for P05.GFP.2 : 37 reads
## Cutoff for P05.MIG.1 : 16 reads
## Cutoff for P05.MIG.2 : 31 reads
## Cutoff for P05.CTR.1 : 18 reads
table(FB.seur$HTO_classification.global)
## 
##  Doublet Negative  Singlet 
##    18932     2130    26824
table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]
## 
##   Doublet  Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 
##     18932      2130      2524      2779      2564      2645      2363      2739 
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 
##      2807      2851      2810      2742
# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.998)
## Cutoff for P28.GFP.1 : 43 reads
## Cutoff for P28.GFP.2 : 48 reads
## Cutoff for P28.MIG.1 : 54 reads
## Cutoff for P28.MIG.2 : 32 reads
## Cutoff for P28.CTR.1 : 15 reads
## Cutoff for P05.GFP.1 : 20 reads
## Cutoff for P05.GFP.2 : 41 reads
## Cutoff for P05.MIG.1 : 18 reads
## Cutoff for P05.MIG.2 : 35 reads
## Cutoff for P05.CTR.1 : 20 reads
table(FB.seur$HTO_classification.global)
## 
##  Doublet Negative  Singlet 
##    17761     2819    27306
table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]
## 
##   Doublet  Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 
##     17761      2819      2629      2865      2495      2643      2284      2839 
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 
##      2849      2943      2909      2850
# Demultiplex cells based on HTO enrichment
FB.seur <- HTODemux(FB.seur, assay = 'HTO', positive.quantile = 0.99)
## Cutoff for P28.GFP.1 : 33 reads
## Cutoff for P28.GFP.2 : 40 reads
## Cutoff for P28.MIG.1 : 41 reads
## Cutoff for P28.MIG.2 : 24 reads
## Cutoff for P28.CTR.1 : 10 reads
## Cutoff for P05.GFP.1 : 15 reads
## Cutoff for P05.GFP.2 : 31 reads
## Cutoff for P05.MIG.1 : 14 reads
## Cutoff for P05.MIG.2 : 26 reads
## Cutoff for P05.CTR.1 : 15 reads
table(FB.seur$HTO_classification.global)
## 
##  Doublet Negative  Singlet 
##    20831     1488    25567
table(FB.seur$hash.ID)[c("Doublet","Negative",rownames(FB))]
## 
##   Doublet  Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 
##     20831      1488      2348      2655      2514      2512      2384      2570 
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 
##      2683      2691      2626      2584
# tags q0.99
#cutoff.FB <- c(33,40,41,24,10,
#               15,31,14,26,15)

# tags1/2/6-10 q0.996; tags3-5 q0.99
#cutoff.FB <- c(38,45,41,24,10,
#               18,37,16,31,18)

# tags1/2/6-10 q0.998; tags3-5 q0.99
cutoff.FB <- c(43,48,41,24,10,
               20,41,18,35,20)

names(cutoff.FB) <- rownames(FB.seur)
cutoff.FB
## P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1 
##        43        48        41        24        10        20        41        18 
## P05.MIG.2 P05.CTR.1 
##        35        20
par(mfrow=c(2,5))

#plot(sort(rowSums(t(FB))),
#     xlab="reordered cells", 
#     ylab="UMI counts", log="xy",
#     main="sum")


plot(sort(t(FB)[,1]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[1])
abline(h=cutoff.FB[1],col = color.FB[1])
plot(sort(t(FB)[,2]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[2])
abline(h=cutoff.FB[2],col = color.FB[2])
plot(sort(t(FB)[,3]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[3])
abline(h=cutoff.FB[3],col = color.FB[3])
plot(sort(t(FB)[,4]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[4])
abline(h=cutoff.FB[4],col = color.FB[4])
plot(sort(t(FB)[,5]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[5])
abline(h=cutoff.FB[5],col = color.FB[5])

plot(sort(t(FB)[,6]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[6])
abline(h=cutoff.FB[6],col = color.FB[6])
plot(sort(t(FB)[,7]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[7])
abline(h=cutoff.FB[7],col = color.FB[7])
plot(sort(t(FB)[,8]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[8])
abline(h=cutoff.FB[8],col = color.FB[8])
plot(sort(t(FB)[,9]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[9])
abline(h=cutoff.FB[9],col = color.FB[9])
plot(sort(t(FB)[,10]),
     xlab="reordered cells", 
     ylab="UMI counts", log="xy",
     main=rownames(FB)[10])
abline(h=cutoff.FB[10],col = color.FB[10])

custom.Demux <- function(FBobj,FBcutoff){
  # define decision matrix
  dd <- FBobj@assays[['HTO']]@counts
  dd <- apply(dd, 2, function(x){
    x > FBcutoff
  })
  x1 <- colSums(dd)
  x1[x1>1] <- "Doublet"
  x1[x1==0] <- "Negative"
  x1[x1==1] <- "Singlet"

  x2 <- x1
  
  bc.slt  <- names(x1)[x1=="Singlet"]
  
  for(bc in bc.slt){
    x2[bc] <- rownames(dd)[dd[,bc]]
  }
  
  FBdf <- data.frame(HTO_classification.global=factor(x1,
                                                      levels = c("Doublet","Negative","Singlet")),
                     hash.ID=factor(x2,
                                    levels = c("Doublet","Negative",rownames(dd))))
  return(FBdf)
}

df.FB <- custom.Demux(FB.seur,cutoff.FB)
dim(df.FB)
## [1] 47886     2
df.FB[1:10,]
##                    HTO_classification.global   hash.ID
## AAACCCAAGAAATCCA-1                   Singlet P05.GFP.1
## AAACCCAAGAATCGTA-1                   Singlet P28.GFP.1
## AAACCCAAGAGCCTGA-1                   Doublet   Doublet
## AAACCCAAGATGCTAA-1                   Doublet   Doublet
## AAACCCAAGCACCTGC-1                   Doublet   Doublet
## AAACCCAAGCGATCGA-1                   Singlet P05.MIG.2
## AAACCCAAGGTAGCCA-1                   Singlet P28.GFP.1
## AAACCCAAGGTTGAGC-1                   Singlet P28.GFP.1
## AAACCCAAGTAAACGT-1                   Doublet   Doublet
## AAACCCAAGTAAGCAT-1                   Doublet   Doublet
table(df.FB)
##                          hash.ID
## HTO_classification.global Doublet Negative P28.GFP.1 P28.GFP.2 P28.MIG.1
##                  Doublet    18820        0         0         0         0
##                  Negative       0     1966         0         0         0
##                  Singlet        0        0      2496      2793      2698
##                          hash.ID
## HTO_classification.global P28.MIG.2 P28.CTR.1 P05.GFP.1 P05.GFP.2 P05.MIG.1
##                  Doublet          0         0         0         0         0
##                  Negative         0         0         0         0         0
##                  Singlet       2732      2594      2721      2725      2836
##                          hash.ID
## HTO_classification.global P05.MIG.2 P05.CTR.1
##                  Doublet          0         0
##                  Negative         0         0
##                  Singlet       2785      2720
FB.seur@meta.data[,c("HTO_classification.global","hash.ID")] <- df.FB
FB.seur@meta.data[1:4,]
##                     orig.ident nCount_RNA nFeature_RNA nCount_HTO nFeature_HTO
## AAACCCAAGAAATCCA-1 plus_SIMPLE        221           10        221           10
## AAACCCAAGAATCGTA-1 plus_SIMPLE         80            8         80            8
## AAACCCAAGAGCCTGA-1 plus_SIMPLE        670           10        670           10
## AAACCCAAGATGCTAA-1 plus_SIMPLE        142            9        142            9
##                    HTO_maxID HTO_secondID HTO_margin  HTO_classification
## AAACCCAAGAAATCCA-1 P05.GFP.1    P05.GFP.2  2.0828322           P05.GFP.1
## AAACCCAAGAATCGTA-1 P28.GFP.1    P28.MIG.1  1.0297322           P28.GFP.1
## AAACCCAAGAGCCTGA-1 P28.GFP.1    P05.MIG.2  0.1030854 P05.MIG.2_P28.GFP.1
## AAACCCAAGATGCTAA-1 P05.CTR.1    P28.MIG.2  0.1557043 P05.CTR.1_P28.MIG.2
##                    HTO_classification.global   hash.ID
## AAACCCAAGAAATCCA-1                   Singlet P05.GFP.1
## AAACCCAAGAATCGTA-1                   Singlet P28.GFP.1
## AAACCCAAGAGCCTGA-1                   Doublet   Doublet
## AAACCCAAGATGCTAA-1                   Doublet   Doublet
## if all using the same q, run this one
#FB.seur$HTO_classification.global <- factor(as.character(FB.seur$HTO_classification.global),
#                                            levels = c("Doublet","Negative","Singlet"))
#FB.seur$hash.ID <- factor(as.character(FB.seur$hash.ID),
#                          levels = c("Doublet","Negative",rownames(FB.seur)))
sl_stat <- table(FB.seur$HTO_classification.global)
barplot(sl_stat,ylim = c(0,32500),col = c("#FF6C91","lightgrey","#7CAE00"),main = "Feature Barcode statistics")
text(x=1:3*1.2-0.5,y=sl_stat+1680,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

RidgePlot

with ridge plots, raw
FB.seur@active.ident <- factor(FB.seur$HTO_maxID,levels=rownames(FB.seur))

RidgePlot(FB.seur, assay = 'HTO', features = rownames(FB.seur[['HTO']]),same.y.lims= TRUE, ncol=5,
          cols = color.FB) 

with ridge plots, filtered
#FB.seur@meta.data$hash.ID.sort <- factor(as.character(FB.seur@meta.data$hash.ID),levels = c("Doublet","Negative",levels(FB.seur@active.ident)))
RidgePlot(FB.seur, assay = 'HTO', features = rownames(FB.seur[['HTO']]), ncol=5,same.y.lims= TRUE,
          cols = c("#FF6C91","lightgrey",color.FB),group.by = "hash.ID") 

tSNE for HTOs

# first, remove negative cells from th object
#FB.seur.subset <- FB.seur

# Calculate a tSNE embedding of the HTO data
#DefaultAssay(FB.seur.subset) <- "HTO"

#FB.seur.subset <- ScaleData(FB.seur.subset, features = rownames(FB.seur.subset), verbose = FALSE)
#FB.seur.subset <- RunPCA(FB.seur.subset, features = rownames(FB.seur.subset), approx = FALSE)
#FB.seur.subset <- RunTSNE(FB.seur.subset, dims = 1:10, perplexity = 500, check_duplicates = FALSE)

#saveRDS(FB.seur.subset, "FB0811.seur.subset.rds")
FB.seur.subset <- readRDS("FB0811.seur.subset.rds")
multiplot(
DimPlot(FB.seur.subset, order = rev(rownames(FB.seur)),
        cols = color.FB) + labs(title = "FB Max_ID"), 
DimPlot(FB.seur.subset, order = rev(c("Negative",rownames(FB.seur),"Doublet")), 
        group.by = 'hash.ID', label = F,
        cols = c("lightgrey",color.FB,"#FF6C91")) + 
  labs(title = "FB Filtered_ID") + theme(plot.title = element_text(hjust = 0)) ,
cols = 1)

stat

FB.seur$HTO_classification.global <- factor(as.character(FB.seur$HTO_classification.global),
                                            levels = c("Doublet","Negative","Singlet"))
Idents(FB.seur) <- "HTO_classification.global"
VlnPlot(FB.seur, features = "nCount_RNA", pt.size = 0, log = TRUE, cols = c("#F8766D","lightgrey","#00BA38"))

VlnPlot(FB.seur, features = "nCount_RNA", pt.size = 0, log = TRUE, cols = c("#F8766D","lightgrey","#00BA38")) + geom_jitter(alpha=0.015)

table(FB.seur@meta.data$hash.ID)
## 
##   Doublet  Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 
##     18820      1966      2496      2793      2698      2732      2594      2721 
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 
##      2725      2836      2785      2720
par(mar=c(6,4,2,3))
sl_stat <- table(FB.seur$hash.ID)
barplot(sl_stat,ylim = c(0,23800),
        col = c("#FF6C91","lightgrey",color.FB),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(FB.seur@meta.data$hash.ID.sort), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+1575,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

GEX

mainly follow https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

GEX.seur <- CreateSeuratObject(counts = GEX,
                               min.cells = 3,
                               min.features = 200,
                               project = "plus_SIMPLE")
GEX.seur
## An object of class Seurat 
## 18371 features across 47886 samples within 1 assay 
## Active assay: RNA (18371 features, 0 variable features)

filtering

MT genes

grep("^mt-",rownames(GEX),value = T)
##  [1] "mt-Nd1"  "mt-Nd2"  "mt-Co1"  "mt-Co2"  "mt-Atp8" "mt-Atp6" "mt-Co3" 
##  [8] "mt-Nd3"  "mt-Nd4l" "mt-Nd4"  "mt-Nd5"  "mt-Nd6"  "mt-Cytb"
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)
GEX.seur[["percent.mt"]] <- PercentageFeatureSet(GEX.seur, features = MT_gene)

RB_gene <- grep("^Rps|^Rpl",rownames(GEX.seur),value = T)
GEX.seur[["percent.rb"]] <- PercentageFeatureSet(GEX.seur, features = RB_gene)

# Visualize QC metrics as a violin plot
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.01)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0)

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb")
plota + plotb + plotc

add FB info

GEX.seur$FB.info <- FB.seur@meta.data[rownames(GEX.seur@meta.data),"hash.ID"]
#head(GEX.seur@meta.data)
table(GEX.seur$FB.info)
## 
##   Doublet  Negative P28.GFP.1 P28.GFP.2 P28.MIG.1 P28.MIG.2 P28.CTR.1 P05.GFP.1 
##     18820      1966      2496      2793      2698      2732      2594      2721 
## P05.GFP.2 P05.MIG.1 P05.MIG.2 P05.CTR.1 
##      2725      2836      2785      2720
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = c("#FF6C91","lightgrey",color.FB),
        pt.size = 0.1, split.by = "FB.info")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = c("#FF6C91","lightgrey",color.FB),
        pt.size = 0, split.by = "FB.info")

GEX.seur <- subset(GEX.seur, subset = FB.info!="Doublet" & FB.info!="Negative")
GEX.seur
## An object of class Seurat 
## 18371 features across 27100 samples within 1 assay 
## Active assay: RNA (18371 features, 0 variable features)
VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.1)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = c("#FF6C91","lightgrey",color.FB),
        pt.size = 0.1, split.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = c("#FF6C91","lightgrey",color.FB),
        pt.size = 0, split.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb") 
plota + plotb + plotc

filtering

GEX.seur <- subset(GEX.seur, subset = percent.mt < 10 & nFeature_RNA > 800 & nFeature_RNA < 3200 & nCount_RNA < 10000)
GEX.seur
## An object of class Seurat 
## 18371 features across 25338 samples within 1 assay 
## Active assay: RNA (18371 features, 0 variable features)

filtered obj

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0.1)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt","percent.rb"), ncol = 4, pt.size = 0)

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FBraw,
        pt.size = 0.1, split.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FBraw,
        pt.size = 0, split.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb") 
plota + plotb + plotc

par(mar=c(6,4,2,3))
sl_stat <- table(FB.seur$hash.ID)
barplot(sl_stat,ylim = c(0,23800),
        col = c("#FF6C91","lightgrey",color.FB),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(FB.seur@meta.data$hash.ID), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+1575,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.info)
barplot(sl_stat,ylim = c(0,4800),
        col = c("#FF6C91","lightgrey",color.FB),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(FB.seur@meta.data$hash.ID), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+345,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

check cellcycle

s.genes <- Hmisc::capitalize(tolower(cc.genes.updated.2019$s.genes))
g2m.genes <- Hmisc::capitalize(tolower(cc.genes.updated.2019$g2m.genes))

GEX.seur <- CellCycleScoring(GEX.seur,
                             s.features = s.genes,
                             g2m.features = g2m.genes)
VlnPlot(GEX.seur, features = c("S.Score", "G2M.Score"), group.by = "FB.info", 
    ncol = 2, pt.size = 0.1)

GEX.seur.cc <- GEX.seur

GEX.seur.cc <- NormalizeData(GEX.seur.cc)
GEX.seur.cc <- FindVariableFeatures(GEX.seur.cc)
GEX.seur.cc <- ScaleData(GEX.seur.cc)
Idents(GEX.seur.cc) <- "Phase"
RidgePlot(GEX.seur.cc, features = c("Pcna", "Top2a", "Mcm6", "Mki67"), ncol = 2)

GEX.seur.cc <- RunPCA(GEX.seur.cc, features = c(s.genes, g2m.genes))
DimPlot(GEX.seur.cc, reduction = 'pca')

seems there’re partly cycling !

Markers and Clusters

Normalizing

GEX.seur <- NormalizeData(GEX.seur, normalization.method = "LogNormalize", scale.factor = 10000)
GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1500)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 200)
##   [1] "Cd74"          "Spp1"          "H2-Aa"         "H2-Eb1"       
##   [5] "H2-Ab1"        "Ccl4"          "Cxcl10"        "Ccl5"         
##   [9] "Pf4"           "Ccl3"          "Spp1-EGFP"     "S100a6"       
##  [13] "Ccl12"         "Hist1h1b"      "Crip1"         "S100a4"       
##  [17] "Mki67"         "Cxcl9"         "Ccl2"          "Cd209a"       
##  [21] "Ifi27l2a"      "Cxcl13"        "Rsad2"         "Top2a"        
##  [25] "Apoe"          "Fn1"           "Xist"          "Lgals3"       
##  [29] "Ccr2"          "Fabp5"         "Plac8"         "Ifitm3"       
##  [33] "Cst7"          "Lpl"           "Lyve1"         "Hist1h2ae"    
##  [37] "Lyz2"          "Cenpf"         "S100a11"       "Mgl2"         
##  [41] "F13a1"         "Hmgb2"         "Hist1h2ap"     "Ms4a7"        
##  [45] "Meg3"          "Snhg11"        "Ch25h"         "Ccl7"         
##  [49] "Ifit3"         "Isg15"         "Cd163"         "Ifi203"       
##  [53] "Prc1"          "Ifit2"         "AA467197"      "Vim"          
##  [57] "Cd209f"        "Stmn1"         "Iigp1"         "Gpnmb"        
##  [61] "Ear2"          "Mrc1"          "Syt1"          "Pclaf"        
##  [65] "Wfdc17"        "Cd52"          "Ifit1"         "Napsa"        
##  [69] "Ube2c"         "Egr1"          "Lgals1"        "Hmmr"         
##  [73] "Mmp12"         "S100a10"       "Plbd1"         "Ahnak"        
##  [77] "Ccl8"          "Nusap1"        "Il1b"          "Apoc1"        
##  [81] "Nkg7"          "Cenpa"         "Cd72"          "Arg1"         
##  [85] "Cd69"          "Cybb"          "Csf1"          "Hmox1"        
##  [89] "Gria2"         "Ccnd2"         "Nrxn3"         "Cd3d"         
##  [93] "Atp1b1"        "Aspm"          "Fxyd5"         "Cenpe"        
##  [97] "Ccl6"          "Ifitm6"        "Birc5"         "Apoc2"        
## [101] "Gm26737"       "Rian"          "Cd9"           "Serpina3g"    
## [105] "Epcam"         "Ctsd"          "Apod"          "Ccr7"         
## [109] "Mndal"         "Cdca8"         "Cd3g"          "Cadps"        
## [113] "Pbld1"         "H2-K1"         "Ptprcap"       "Igf1"         
## [117] "H2afz"         "Ifi204"        "Ms4a4b"        "Hist1h3c"     
## [121] "Tpx2"          "Spn"           "Ptgds"         "Olfm1"        
## [125] "Ndn"           "Lgals3bp"      "Slamf7"        "Trbc2"        
## [129] "Gdf15"         "Il12b"         "Fth1"          "H2-D1"        
## [133] "Mif"           "Bst2"          "Klrd1"         "Igfbp5"       
## [137] "Clec4b1"       "Slpi"          "Stmn3"         "Gad1"         
## [141] "Smc2"          "Ly6a"          "Saa3"          "Ctsb"         
## [145] "Clec10a"       "Nfasc"         "Ptprn"         "Kif11"        
## [149] "Kctd16"        "Slfn1"         "B230312C02Rik" "Ifitm2"       
## [153] "Stmn2"         "Acod1"         "Ttr"           "Gzmb"         
## [157] "Kctd14"        "Esco2"         "Ifi205"        "Slfn5"        
## [161] "Igfbp2"        "Hp"            "Cd3e"          "H2afx"        
## [165] "Rcan1"         "Cdkn1a"        "Knl1"          "Smc4"         
## [169] "Chl1"          "Ifi211"        "Ace"           "Il2rb"        
## [173] "Gimap3"        "Ccnb1"         "Rrm2"          "Ftl1"         
## [177] "Hist1h1e"      "Ank3"          "Cd24a"         "Traf1"        
## [181] "Adarb2"        "B830012L14Rik" "Pbk"           "Emb"          
## [185] "Tmsb10"        "Gad2"          "Hba-a2"        "Dab2"         
## [189] "Cd209g"        "B2m"           "Gimap7"        "Hba-a1"       
## [193] "Ms4a4c"        "Trac"          "Cldn11"        "Cd36"         
## [197] "Cd63"          "Jaml"          "Runx3"         "Cd8b1"
GEX.seur <- ScaleData(GEX.seur, features = rownames(GEX.seur))
## Centering and scaling data matrix

PCA

# exclude MT genes  and more 

#DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^Gm|^Hist",rownames(GEX.seur),value = T)

DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AY|^AC|^AA|^AW|^BC|^Gm|^Hist|Rik$|-ps",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))


VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))
## PC_ 1 
## Positive:  Cst7, Ifitm3, Ifi27l2a, Ifi204, Isg15, B2m, Cd72, Lgals3bp, Fth1, Ch25h 
##     Cd52, Bst2, H2-K1, Ifi207, Ifit2, Oasl2, Ctsb, Fxyd5, Ifit3, H2-D1 
##     Zbp1, Rsad2, Ccl3, Irf7, Ctsd, Ctsz, Slfn5, Tspo, Ccl2, Cd9 
## Negative:  Crybb1, Ccr5, Sox4, Fchsd2, Mrc1, Igfbp4, Rgs7bp, Hmgb1, Cryba4, Cbfa2t3 
##     Il7r, Gbp7, Gas6, Atrx, Stab1, Ccr1, Zbtb20, Khdrbs3, Tmem52, Upk1b 
##     Slc12a2, Stox2, Bod1l, Lsp1, H1f0, Clec4a1, Smc3, Lcorl, Ednrb, Cask 
## PC_ 2 
## Positive:  Prc1, Pclaf, Pbk, Knl1, Kif15, Esco2, Aspm, Racgap1, Spc24, Lockd 
##     Mis18bp1, Stmn1, Smc2, Bub1b, Neil3, Ccna2, Sgo2a, Ccnb1, Ccnf, Plk1 
##     H2afx, Spc25, Cit, Tk1, Ncapg, Shcbp1, Diaph3, Kif4, Trim59, Kif22 
## Negative:  Ctsd, Tyrobp, Ctsb, Timp2, Apoe, Ctsl, Cd63, B2m, Abca1, H2-D1 
##     H2-K1, Cd9, Ctsz, Mpeg1, Fth1, Lgals3bp, Grn, Npc2, Ftl1, Lag3 
##     Lyz2, Pmp22, Xist, H2-T23, Slfn5, Ly86, Cst7, Ifi204, Ccl9, Abcg1 
## PC_ 3 
## Positive:  Ctsd, Cd9, Ctsl, Ctsb, Ctsz, Syngr1, Cst7, Cd63, Serpine2, Lpl 
##     Tyrobp, Grn, Timp2, Ly86, Ccl3, Fabp5, Mif, Baiap2l2, Pld3, Atp6v0c 
##     Hif1a, Csf1, Spp1-EGFP, Pmp22, Gnas, Cadm1, Myo1e, Spp1, Fam20c, Cpd 
## Negative:  H2-Aa, H2-Eb1, S100a6, Ifi203, Plbd1, H2-Ab1, S100a11, Cd74, S100a4, Mndal 
##     Crip1, Ccr2, Cybb, Vim, Napsa, Ahnak, Clec12a, Slamf7, Itga4, S100a10 
##     Emb, Mgl2, Ms4a4c, Clec10a, Ciita, Ifi211, Wfdc17, Ifitm2, Ifit3b, Ifit3 
## PC_ 4 
## Positive:  Lgals3, Crip1, Lgals1, Vim, Fxyd5, H2-Ab1, Ftl1, H2-Eb1, Anxa2, H2-Aa 
##     S100a6, S100a11, S100a4, Cd74, Igf1, Plbd1, Iqgap1, Ifitm2, Ccr2, Anxa5 
##     Alcam, Ahnak, Wfdc17, Napsa, Lpl, Itga4, Fabp5, Bhlhe40, Emp3, Clec12a 
## Negative:  Ifit3, Ifit2, Ifit1, Rsad2, Iigp1, Ifit3b, Isg15, Ifi204, Cxcl10, Slfn5 
##     Usp18, Parp14, Oasl2, Phf11d, Ifi206, Ifi209, Herc6, Ccl12, Stat1, Ifi213 
##     Oasl1, Fgl2, Rtp4, Ifi208, Slfn8, Ifi211, Trim30a, Irf7, Gbp2, Tor3a 
## PC_ 5 
## Positive:  H2-Ab1, H2-Aa, H2-Eb1, S100a11, S100a6, Cd74, Klrd1, Plbd1, Ccr2, S100a4 
##     Cd9, Napsa, Vim, Crip1, Slamf7, S100a10, H2-DMa, Il7r, Mcemp1, Csmd3 
##     Alcam, Ciita, Ctsd, Cytip, H2-DMb1, Emb, Olfm1, Anxa2, Cd34, Cd209a 
## Negative:  Pf4, Ms4a7, Mrc1, F13a1, Cbr2, Dab2, Pla2g7, Cd163, Lyve1, Apoe 
##     Cp, Ednrb, Gas6, Folr2, Clec4n, Colec12, Igfbp4, Cd36, Pltp, Stab1 
##     Igf1, Cd38, Itsn1, Msr1, Clec4a1, Gpx3, Myo5a, Clec10a, Adam33, Mgl2
length(VariableFeatures(GEX.seur))
## [1] 1229
head(VariableFeatures(GEX.seur),300)
##   [1] "Cd74"      "Spp1"      "H2-Aa"     "H2-Eb1"    "H2-Ab1"    "Ccl4"     
##   [7] "Cxcl10"    "Ccl5"      "Pf4"       "Ccl3"      "Spp1-EGFP" "S100a6"   
##  [13] "Ccl12"     "Crip1"     "S100a4"    "Cxcl9"     "Ccl2"      "Cd209a"   
##  [19] "Ifi27l2a"  "Cxcl13"    "Rsad2"     "Apoe"      "Fn1"       "Xist"     
##  [25] "Lgals3"    "Ccr2"      "Fabp5"     "Plac8"     "Ifitm3"    "Cst7"     
##  [31] "Lpl"       "Lyve1"     "Lyz2"      "S100a11"   "Mgl2"      "F13a1"    
##  [37] "Ms4a7"     "Meg3"      "Snhg11"    "Ch25h"     "Ccl7"      "Ifit3"    
##  [43] "Isg15"     "Cd163"     "Ifi203"    "Prc1"      "Ifit2"     "Vim"      
##  [49] "Cd209f"    "Stmn1"     "Iigp1"     "Gpnmb"     "Ear2"      "Mrc1"     
##  [55] "Syt1"      "Pclaf"     "Wfdc17"    "Cd52"      "Ifit1"     "Napsa"    
##  [61] "Egr1"      "Lgals1"    "Mmp12"     "S100a10"   "Plbd1"     "Ahnak"    
##  [67] "Ccl8"      "Il1b"      "Apoc1"     "Nkg7"      "Cd72"      "Arg1"     
##  [73] "Cd69"      "Cybb"      "Csf1"      "Hmox1"     "Gria2"     "Ccnd2"    
##  [79] "Nrxn3"     "Cd3d"      "Atp1b1"    "Aspm"      "Fxyd5"     "Ccl6"     
##  [85] "Ifitm6"    "Apoc2"     "Rian"      "Cd9"       "Serpina3g" "Epcam"    
##  [91] "Ctsd"      "Apod"      "Ccr7"      "Mndal"     "Cd3g"      "Cadps"    
##  [97] "Pbld1"     "H2-K1"     "Ptprcap"   "Igf1"      "H2afz"     "Ifi204"   
## [103] "Ms4a4b"    "Spn"       "Ptgds"     "Olfm1"     "Ndn"       "Lgals3bp" 
## [109] "Slamf7"    "Gdf15"     "Il12b"     "Fth1"      "H2-D1"     "Mif"      
## [115] "Bst2"      "Klrd1"     "Igfbp5"    "Clec4b1"   "Slpi"      "Stmn3"    
## [121] "Gad1"      "Smc2"      "Ly6a"      "Saa3"      "Ctsb"      "Clec10a"  
## [127] "Nfasc"     "Ptprn"     "Kctd16"    "Slfn1"     "Ifitm2"    "Stmn2"    
## [133] "Acod1"     "Ttr"       "Gzmb"      "Kctd14"    "Esco2"     "Ifi205"   
## [139] "Slfn5"     "Igfbp2"    "Hp"        "Cd3e"      "H2afx"     "Rcan1"    
## [145] "Cdkn1a"    "Knl1"      "Chl1"      "Ifi211"    "Ace"       "Il2rb"    
## [151] "Gimap3"    "Ccnb1"     "Ftl1"      "Ank3"      "Cd24a"     "Adarb2"   
## [157] "Pbk"       "Emb"       "Tmsb10"    "Gad2"      "Dab2"      "Cd209g"   
## [163] "B2m"       "Gimap7"    "Ms4a4c"    "Cldn11"    "Cd36"      "Cd63"     
## [169] "Jaml"      "Runx3"     "Cd8b1"     "Kit"       "Ciita"     "Ncam1"    
## [175] "Adgre5"    "Ifi207"    "Grik1"     "Cd38"      "Lilrb4a"   "Itgal"    
## [181] "Cd5l"      "Nr4a3"     "Clstn3"    "Kif15"     "Plp1"      "Mis18bp1" 
## [187] "Ly6h"      "Ctsz"      "Slc6a1"    "Krt80"     "Gbp2"      "Grm5"     
## [193] "Lsp1"      "Clec12a"   "Upb1"      "Enpp2"     "Dqx1"      "Ifi44"    
## [199] "Mt1"       "Ccl9"      "Syngr1"    "Cbr2"      "Sh2d2a"    "Tnr"      
## [205] "H2-Q7"     "Itga4"     "Snca"      "Acp5"      "Id2"       "Cd226"    
## [211] "Drd4"      "Mcemp1"    "Ifit3b"    "Ccr1"      "Mcf2l"     "Nsg1"     
## [217] "Pla2g2d"   "Tspo"      "Cxcl2"     "Adgre4"    "Ccl24"     "Irf7"     
## [223] "Marco"     "C4b"       "Cspg4"     "Bex2"      "Cd93"      "Bub1b"    
## [229] "Map1b"     "Racgap1"   "Nap1l5"    "Sgo2a"     "Nap1l3"    "Dennd5b"  
## [235] "Wdcp"      "Ctsl"      "Adamdec1"  "Gapdh"     "Opcml"     "Atf3"     
## [241] "Tc2n"      "Oasl2"     "Anxa1"     "Ctla2a"    "Cmpk2"     "Gnas"     
## [247] "Klrb1c"    "Rbms3"     "Pcdh17"    "Aldh1a3"   "Scn2a"     "Vstm2l"   
## [253] "Cstb"      "Adcy1"     "Clec4n"    "Sst"       "Spc24"     "Peg3"     
## [259] "Gap43"     "Baiap2l2"  "Brinp2"    "Ppp1r14a"  "Rnase2a"   "Ctsw"     
## [265] "Sct"       "Ifitm1"    "Fpr1"      "Ms4a6c"    "Phf11b"    "Ank"      
## [271] "Nfib"      "Kmo"       "Myh4"      "Lmnb1"     "Tnfsf18"   "Cox6a2"   
## [277] "Epha5"     "Kcnq1ot1"  "Emp3"      "Fabp4"     "Cd7"       "Kcnn2"    
## [283] "Sh3tc2"    "Gfod2"     "Capg"      "Cadm2"     "Tk1"       "Fpr2"     
## [289] "Tsix"      "Auts2"     "Lockd"     "Spc25"     "Ly86"      "Ntm"      
## [295] "Klrb1b"    "Amph"      "Cxcl16"    "Blk"       "Ifit1bl1"  "Timp2"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB)

DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 25338
## Number of edges: 669263
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7603
## Number of communities: 20
## Elapsed time: 5 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 50, seed.use = 287)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 11:40:13 UMAP embedding parameters a = 0.9922 b = 1.112
## 11:40:13 Read 25338 rows and found 20 numeric columns
## 11:40:13 Using Annoy for neighbor search, n_neighbors = 50
## 11:40:13 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 11:40:16 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\Rtmpe8iaop\filee42c639d6904
## 11:40:16 Searching Annoy index using 1 thread, search_k = 5000
## 11:40:29 Annoy recall = 100%
## 11:40:30 Commencing smooth kNN distance calibration using 1 thread
## 11:40:33 Initializing from normalized Laplacian + noise
## 11:40:34 Commencing optimization for 200 epochs, with 1861048 positive edges
## 11:41:08 Optimization finished
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'tsne', group.by = "FB.info", split.by = "FB.info", 
        ncol = 5, cols = color.FB)

DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.info", split.by = "FB.info", 
        ncol = 5, cols = color.FB)

GEX.seur$cnt <- factor(gsub(".1$|.2$","",as.character(GEX.seur$FB.info)),
                       levels = c("P05.CTR","P05.MIG","P05.GFP",
                                  "P28.CTR","P28.MIG","P28.GFP"))
color.cnt <- color.FB[c(10,9,7,
                        5,4,1)]
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)

check some markers

markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
                 "Cd74","Lyz2","Ccl4",
                 "Aif1","P2ry12","C1qa","Spp1",
                 "Top2a","Pcna","Mki67","Mcm6",
                 "Cx3cr1","Il4ra","Il13ra1","Spp1-EGFP",
                 "Fabp5","Hmox1","Ms4a7","Cenpa")
FeaturePlot(GEX.seur, 
            features = markers.mig,
            ncol = 4)

DotPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
                               "Fcer1g","Il4ra","Il13ra1"))

VlnPlot(GEX.seur, features = c("Cx3cr1","Spp1","Aif1", # Aif1: Iba1
                               "Fcer1g","Il4ra","Il13ra1"))

GEX.seur@meta.data$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                           levels = c(1,2,7,5,
                                                      8,14,13,16,
                                                      9,
                                                      0,4,6,10,3,17,
                                                      11,12,
                                                      
                                                      15,18,19))
DotPlot(GEX.seur, features = markers.mig, group.by = "sort_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.info,
      clusters=GEX.seur$sort_clusters)[c(3:12),],
                   main = "Cell Count",
                   gaps_row = c(2,4,5,7,9),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.info,
                                clusters=GEX.seur$sort_clusters)[c(3:12),]),
                   main = "Cell Ratio",
                   gaps_row = c(2,4,5,7,9),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"

GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, 
                                  min.pct = 0.05,
                                  test.use = "MAST",
                                  logfc.threshold = 0.25)
#GEX.markers.pre <- read.table("sc10x_LYN.marker.sort0811.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## # A tibble: 160 x 7
## # Groups:   cluster [20]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene         
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>        
##  1 6.21e-297      0.422 0.995 0.988 1.14e-292 1       Rps11        
##  2 1.47e-280      0.533 0.961 0.89  2.70e-276 1       Ltc4s        
##  3 1.17e-199      0.453 0.966 0.928 2.15e-195 1       Maf          
##  4 8.30e-180      0.548 0.816 0.677 1.52e-175 1       Crybb1       
##  5 2.76e-174      0.529 0.817 0.722 5.07e-170 1       Ccr5         
##  6 2.12e-146      0.502 0.654 0.528 3.89e-142 1       2410006H16Rik
##  7 1.84e-108      0.426 0.648 0.56  3.38e-104 1       Fchsd2       
##  8 8.56e-105      0.426 0.428 0.33  1.57e-100 1       Fcgrt        
##  9 0              0.812 0.91  0.669 0         2       Crybb1       
## 10 1.30e-206      0.672 0.454 0.199 2.39e-202 2       Igfbp4       
## # ... with 150 more rows
GEX.markers.sort <- GEX.markers.pre
GEX.markers.sort$cluster <- factor(as.character(GEX.markers.sort$cluster),
                          levels = levels(GEX.seur$sort_clusters))


markers.pre_t48 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.1) %>%
                   top_n(n = 48, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t60 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[1:64]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[65:128]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[129:192]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[193:256]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[257:320]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[321:384]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[385:448]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[449:512]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[513:576]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[577:640]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[641:714]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

check.ref1 <- c("Tmem119","Hexb","Slc2a5","P2ry12",
                "Siglech","Trem2",  # microglia
                "Mrc1","Lyve1","Cd163","Siglec1",  # CAM
                "Ly6c2","Ccr2","Plac8","Anxa8",
                "Nr4a1",  # Monocyte-derived
                "Flt3","Zbtb46","Batf3","Itgae",
                "Clec9a",  # 
                "Tubb3","Actl6b" # locally addede neuron markers
                )

DotPlot(GEX.seur, features = check.ref1, group.by = "sort_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
## Warning in FetchData(object = object, vars = features, cells = cells): The
## following requested variables were not found: Anxa8

C15: BAM (Brain Associated Macrophages)
C18: DC
C19: Neuron

others: Microglia

c("Xist","Tsix") %in% VariableFeatures(GEX.seur)
## [1] TRUE TRUE
c("Xist","Tsix","Rps4x") %in% VariableFeatures(GEX.seur)
## [1]  TRUE  TRUE FALSE
c("Rps4y1","Eif1ay","Ddx3y") %in% VariableFeatures(GEX.seur)
## [1] FALSE FALSE FALSE
c("Xist","Tsix","Uty","Ddx3y","Eif2s3y","Kdm5d") %in% VariableFeatures(GEX.seur)
## [1]  TRUE  TRUE FALSE FALSE FALSE FALSE
markers.Sci2019 <- c("Tmem119","Hexb","Slc2a5","P2ry12","Siglech","Trem2","Sparc","Olfml3",  # Microglia
                     "Mrc1","Lyve1","Cd163","Siglec1","Pf4","Ms4a7","Cbr2",  # CAMs
                     "Ly6c2","Ccr2","Plac8","Anxa8","Nr4a1","Mertk","Zbtb26","Cd209a",  # Monocyte-derived
                     "Flt3","Zbtb46","Batf3","Itgae","Clec9a",  # DCs
                     "Tubb3"  # local added
                     )
FeaturePlot(GEX.seur, 
            features = markers.Sci2019,
            ncol = 4)
## Warning in FetchData(object = object, vars = c(dims, "ident", features), : The
## following requested variables were not found: Anxa8

check.Cell2017 <- c("Hexb","Cst3","P2ry12","Cx3cr1",
                    "Cd9","Ctsd","Cst7","Lpl",
                    "Mrc1","Cd163","S100a4","Cd74",
                    "Cd79b","Rag1","Trbc2","Nkg7",
                    "S100a9","Camp"
                    )

DotPlot(GEX.seur, features = check.Cell2017, group.by = "sort_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
## Warning in FetchData(object = object, vars = features, cells = cells): The
## following requested variables were not found: Rag1, S100a9, Camp

DoubletFinder

library(DoubletFinder)
## Loading required package: KernSmooth
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
## Loading required package: ROCR

## NULL
## [1] "Creating 8446 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## Centering and scaling data matrix
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."
## [1] "Creating 8446 artificial doublets..."
## [1] "Creating Seurat object..."
## [1] "Normalizing Seurat object..."
## [1] "Finding variable genes..."
## [1] "Scaling data..."
## [1] "Running PCA..."
## [1] "Calculating PC distance matrix..."
## [1] "Computing pANN..."
## [1] "Classifying doublets.."

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters", split.by = "DoubletFinder0.05")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters", split.by = "DoubletFinder0.1")

table(GEX.seur@meta.data[,c("sort_clusters","DoubletFinder0.05")])
##              DoubletFinder0.05
## sort_clusters Doublet Singlet
##            1       38    4348
##            2      286    3100
##            7       47     949
##            5      119    1613
##            8       65     825
##            14      15     236
##            13      17     385
##            16      12     151
##            9       71     645
##            0      122    4795
##            4      103    1883
##            6       93    1348
##            10      38     557
##            3      150    2004
##            17      12     115
##            11      30     448
##            12      33     382
##            15      12     164
##            18       2     110
##            19       2      13
table(GEX.seur@meta.data[,c("sort_clusters","DoubletFinder0.1")])
##              DoubletFinder0.1
## sort_clusters Doublet Singlet
##            1       67    4319
##            2      450    2936
##            7      101     895
##            5      234    1498
##            8      137     753
##            14      50     201
##            13      94     308
##            16      39     124
##            9      111     605
##            0      204    4713
##            4      152    1834
##            6      143    1298
##            10      83     512
##            3      311    1843
##            17      63      64
##            11      41     437
##            12      75     340
##            15      98      78
##            18      77      35
##            19       4      11

preAnno

# preAnno1
GEX.seur$preAnno1 <- as.character(GEX.seur$seurat_clusters)


GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(1)] <- "MIG.a1"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(2)] <- "MIG.a2"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(7)] <- "MIG.a3"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(9)] <- "MIG.a4"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(5)] <- "MIG.a5"

GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(8)] <- "MIG.a6"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(14)] <- "MIG.a7"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(13)] <- "MIG.a8"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(16)] <- "MIG.a9"

GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(0)] <- "MIG.b1"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(11)] <- "MIG.b2"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(4)] <- "MIG.b3"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(6)] <- "MIG.b4"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(12)] <- "MIG.b5"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(10)] <- "MIG.b6"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(3)] <- "MIG.b7"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(17)] <- "MIG.b8"

GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(18)] <- "DC"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(15)] <- "BAM"
GEX.seur$preAnno1[GEX.seur$preAnno1 %in% c(19)] <- "Neuron"

GEX.seur$preAnno1 <- factor(GEX.seur$preAnno1,
                            levels = c(paste0("MIG.a",1:9),
                                       paste0("MIG.b",1:8),
                                       "DC","BAM","Neuron"))

# preAnno2
GEX.seur$preAnno2 <- as.character(GEX.seur$preAnno1)

GEX.seur$preAnno2 <- gsub("1|2|3|4|5|6|7|8|9","",GEX.seur$preAnno2)
GEX.seur$preAnno2 <- factor(GEX.seur$preAnno2,
                            levels = c("MIG.a","MIG.b",
                                       "DC","BAM","Neuron"))

# preAnno
GEX.seur$preAnno <- as.character(GEX.seur$preAnno2)

GEX.seur$preAnno <- gsub(".a$|.b$","",GEX.seur$preAnno)
GEX.seur$preAnno <- factor(GEX.seur$preAnno,
                           levels = c("MIG","DC","BAM","Neuron"))
pp.raw.1a <- DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "sort_clusters")
pp.raw.1a

pp.raw.1b <- DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt)
pp.raw.1b

pp.raw.1c <- DimPlot(GEX.seur, reduction = "umap", group.by = "cnt", split.by = "cnt", cols = color.cnt, ncol = 3)
pp.raw.1c 

DimPlot(GEX.seur, reduction = "umap", label = T,label.size = 2.4, group.by = "preAnno1") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno2")

#saveRDS(GEX.seur, "GEX0811.seur.preAnno.rds")
cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2),
  DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB),
rel_widths = c(5,5.6),
ncol = 2)

pp.test <- FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 1:2, split.by = "FB.info", combine = F)
cowplot::plot_grid(
  plotlist = pp.test,
  ncol = 5
)

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", split.by = "FB.info", cols = color.FB, ncol = 5)

cowplot::plot_grid(
FeaturePlot(GEX.seur, features = "Spp1", reduction = "pca",dims = 3:4),
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB),
rel_widths = c(5,5.6),
ncol = 2)

cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$cnt,
      clusters=GEX.seur$preAnno),
                   main = "Cell Count",
                   #gaps_row = c(2,4,5,7,9),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$cnt,
                                clusters=GEX.seur$preAnno)),
                   main = "Cell Ratio",
                   #gaps_row = c(2,4,5,7,9),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

table(GEX.seur@meta.data[,c("cnt","preAnno")])
##          preAnno
## cnt        MIG   DC  BAM Neuron
##   P05.CTR 2598    2   12      5
##   P05.MIG 5362   11   32      5
##   P05.GFP 5126    1   32      2
##   P28.CTR 2372    3   14      0
##   P28.MIG 4980   36   27      2
##   P28.GFP 4597   59   59      1
rowRatio(table(GEX.seur@meta.data[,c("cnt","preAnno")]))
##          preAnno
## cnt                MIG           DC          BAM       Neuron
##   P05.CTR 0.9927397784 0.0007642339 0.0045854031 0.0019105846
##   P05.MIG 0.9911275416 0.0020332717 0.0059149723 0.0009242144
##   P05.GFP 0.9932183685 0.0001937609 0.0062003488 0.0003875218
##   P28.CTR 0.9928840519 0.0012557555 0.0058601925 0.0000000000
##   P28.MIG 0.9871159564 0.0071357780 0.0053518335 0.0003964321
##   P28.GFP 0.9747667515 0.0125106022 0.0125106022 0.0002120441
marker.fig0921 <- c("Ptprc","P2ry12","Spp1","Spp1-EGFP",
                    "Cx3cr1","Csf1r","Aif1","C1qa",
                    "Ccl4","Lpl","Fabp5","Csf1",
                    "Ctsf","Ccr5","Arsb","Itgb1",
                    "Cd81","Lag3","Cd34","H2-K1",
                    "Top2a","Pcna","Mki67","Mcm6",
                    "Cd74","H2-Aa","Ccr2","Cd209a",
                    "Mrc1","Pf4","Dab2","Lyve1",
                    "Tubb3","Actl6b","Syt1","Sst"
                    )

pp.raw.2a <- FeaturePlot(GEX.seur, features = marker.fig0921, ncol = 4,raster = T, pt.size = 2.75)
pp.raw.2a

pp.raw.2b <- DotPlot(GEX.seur, features = marker.fig0921, group.by = "preAnno1",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)
pp.raw.2b